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Abstract 

A  multivariate  dynamic  Markov  model  is  formulated  to 
describe  the  possible  effect  of  a  generic  chemical  toxin  on  a  generic 
cell  population  in  an  organ.  Asymptotic  methods  (large  cell 
population)  are  used  to  show  that  numbers  and  toxin  may  be 
jointly  Gaussian/normally  distributed;  the  joint  stochastic  process 
is  approximately  Ornstein-Uhlenbeck.  Application  to  dose- 
response,  and  hence  risk  analysis,  is  briefly  discussed. 


EXECUTIVE  SUMMARY 

This  paper  presents  a  generic  model  for  the  influence  upon  the  cells  in  an 
c:  ui  caused  by  dosage  by  a  chemical  toxin.  The  model  is  presently  non-specific 
to  either  organ  or  toxin,  but  provides  a  starting  point  for  specialization.  The 
model  is  probabilistic/stochastic  without  requiring  computer  simulation  (it  relies 
on  asymptotics  appropriate  for  organs  containing  a  large  number  of  interactive 
cells). 

The  organ  disfunctionality  associated  with  toxin  input  that  is  modeled  is  the 
(probability  distribution  of  the)  number  of  functionally  active  cells  present:  that 
number  tends  to  be  reduced  by  toxin  input.  Inter-cell  signalling  is  modeled  by  a 
pairing  device;  this  functionality  is  due  for  later  elaboration.  The  effects  of  cell 
death  by  necrosis  and  apoptosis  are  modeled.  The  present  model  omits  explicit 
consideration  of  such  issues  as  spatiality  of  cells  within  organs,  enzyme  function 
and  other  important  features. 

Application  of  the  present  model  is  made  to  consideration  of  the  shape  of  a 
dose-response  function  at  low  doses. 


AN  EXPLORATORY  STOCHASTIC  MODEL 

FOR  TOXIC  EFFECTS  ON  CELLS 

R.  L.  Carpenter 

D.  P.  Gaver 

P.  A.  Jacobs 

1.   Introduction 

Application  of  laboratory  toxicology  data  to  environmental  and  human 
problems  of  risk  assessment  almost  always  requires  extrapolation  of  the  data 
from  the  experimentally  used  dose  regimen  to  the  exposure  conditions  of 
practical  concern,  and  from  the  animal  species  tested  to  the  species  of  concern 
(usually  man).  This  extrapolation  and  estimation  process  is  known  as  chemical 
risk  assessment.  The  risk  assessment  process  has  undergone  considerable 
evolution  in  the  last  ten  years,  moving  from  a  qualitative  basis  for  decision 
making  to  an  increasingly  quantitative  basis,  and  from  the  use  of  default 
assumptions  to  the  application  of  mathematical  models  as  tools  upon  which  to 
base  decisions.  In  the  context  of  determining  safe  human  exposure  limits  to 
potentially  toxic  chemicals,  there  are  two  sub-tasks  to  be  accomplished: 
estimation  of  low  dose  risk  in  animals  and  the  subsequent  conversion  of  animal 
risk  estimates  to  human  risk  estimates. 

Dose  extrapolation,  the  first  sub-task,  can  be  accomplished  either  by 
assuming  biological  response  varies  linearly  with  dose,  or  by  using 
physiologically-based  pharmacokinetic  (PBPK)  models.  Multi-compartment 
physiological  models  are  formulated  using  actual  tissue  volumes  from  the 
experimental  and  target  species  and  actual  perfusion  rates  to  provide  for 


chemical  transport  between  the  compartments.  Thus  the  pharmacokinetics  of 
high  to  low-dose  extrapolation  become  amenable  to  calculation,  and  external 
measures  of  dose  can  be  translated  to  concentrations  in  the  target  organ  (internal 
dose).  Target-organ  chemical  concentration  may  be  translated  into  estimates  of 
risk  if  a  mechanistically-based  model  can  be  used  to  relate  chemical  presence  in 
the  organ  to  harmful  outcome.  Such  models  exist  for  carcinogens,  in  the  form  of 
the  widely-used  linearized  multistage  model.  This  model  is  based  on  the 
generalized  concept  that  chemical  alteration  of  the  cellular  genes  may  give  rise  to 
permanent,  heritable  changes  in  the  genetic  information  stored  in  the  cell  nucleus 
and  lead  to  phenotypic  changes  in  the  altered  cells  that  ultimately  cause  the 
formation  of  malignant  tumors. 

Unfortunately,  no  similar  models  exist  for  chemical  toxic  effects  other  than 
carcinogenesis.  The  mathematical  expression  of  even  the  relatively  simple 
concept  of  genetic  alteration  leading  to  cancer  involves  significant  simplification 
of  biological  reality  and  significant  mathematical  complexity.  This  state  of  affairs 
is  exacerbated  when  one  attempts  to  accurately  describe  the  interactions  leading 
to  loss  of  cell  function  and  cell  death  in  tissues  of  a  whole  animal.  The  multi- 
layered  control  and  response  systems  present  in  an  intact  living  animal  are 
poorly  understood  and  have  not  been  modeled  as  yet.  Nonetheless,  these  control 
mechanisms  defend  against  the  majority  of  toxic  effects  observed  as  a  result  of 
chemical  exposure,  and  their  failure  represents  most  of  what  is  observed  as  the 
expression  of  toxicity. 

This  paper  describes  initial  formulations  of  models  describing  cellular 
response  to  toxic  insult.  Our  model  formulation  takes  into  account  aspects  of  the 
current  state  of  knowledge  concerning  the  control  and  regulation  of  cellular 


properties  in  tissue.  The  design  of  this  model  is  not  to  provide  a  detailed 
description  of  the  response  of  specific  organs  to  toxic  insult,  but  rather  to  begin 
the  process  of  mathematically  describing  the  recent  significant  advances  in 
knowledge  that  have  been  made  in  cell  biology,  hormonal  regulation,  and 
hormetic  control  mechanisms.  Considerable  extension  of  model  formulation  will 
be  needed  to  apply  it  to  specific  organs,  since  the  complex  geometry  of  organ 
architecture  is  not  present.  To  have  included  complex,  three-dimensional 
relationships  between  different  cell  types  would  have  increased  the  mathematical 
complexity  considerably.  The  current  formulation  of  the  model  is  suitable  for 
experimental  evaluation  in  cell  culture  if  the  experimental  conditions  are 
properly  chosen. 

2.  Biological  Background 

The  tissue  making  up  organs  almost  always  involves  a  complex  geometrical 
juxtaposition  of  several  cell  types  having  specific  and  often  overlapping 
functions.  This  tissue  architecture  is  maintained  by  the  presence  of  a  nonliving 
matrix  of  proteins  collectively  known  as  a  basement  membrane.  Interactions 
between  this  membrane  and  the  cell  are  known  to  be  critical  to  its  stability  as  a 
mature,  functioning  entity.  The  whole  of  the  tissue  is  permeated  by  branching 
blood  vessels,  each  generation  of  which  is  successively  smaller;  these  serve  to 
provide  a  constant,  stable  milieu  in  which  the  cells  exist.  Nutrients,  control 
signals  in  the  form  of  specific  biomolecules,  and  xenobiotic  chemicals  are  brought 
to  the  immediate  surroundings  of  the  cell  by  the  vasculature  in  the  tissue  and 
cellular  metabolic  products;  the  products  of  energy  metabolism  and  cellularly 
derived  control  biomolecules  are  removed  from  the  cellular  microenvironment 
by  the  same  means. 


The  state  of  the  cell  at  any  time  is  a  reflection  of  its  age,  the  summation  of  the 
control  chemicals  reaching  and  leaving  the  cell,  the  effects  of  xenobiotic 
chemicals  present,  if  any,  and  the  state  of  the  cells  surrounding  it.  Cellular 
contact  with  the  surrounding  cells  and  basement  membrane  act  to  supply 
chemical  signals  that  mod.  late  its  activity.  A  given  cell  may  be  (i)  nearly 
quiescent,  (ii)  active  biochemically,  i.e.,  producing  metabolites  of  absorbed 
materials  for  its  own  use  or  for  export,  (iii)  in  a  state  of  stress  due  to  shortage  or 
excess  of  biochemical  molecules,  (iv)  in  a  process  of  programmed  cell  death 
(apoptosis),  (v)  in  the  process  of  dying  from  chemical  insult  (necrosis),  or  (vi) 
dividing  to  form  progenitor  cells  in  response  to  a  need  to  replace  cells  already 
lost;  these  conditions  need  not  be  mutually  exclusive  at  any  point  in  time.  Some 
specific  cells  may  alter  or  completely  change  their  observable  characteristics  in 
response  to  chemical  signals.  The  most  well  known  of  these  cells  are  the 
pluripotent  stem  cells  of  the  hemopoetic  system. 

3.   Prototypic  Mathematical  Models  for  Cellular  Response  to  a  Chemical 
Toxin 

In  this  section  we  propose  several  simple  models  to  describe  the  interaction 
between  cells  and  a  chemical  toxin  with  which  they  are  in  contact.  The  aim  is  to 
describe  anti-toxin  functionality  in  a  parsimonious  manner,  so  certain  details  are 
omitted  that  can  feature  in  later  models. 

Consider  a  collection  of  interacting  cells  in  proximity  in  an  organ  such  as  the 
liver.  In  particular,  these  may  be  the  cells  that  line  the  sinusoids  in  the  liver,  the 
function  of  which  is  to  remove  waste  from  the  blood  flowing  through. 

Focus  first  on  the  cell  cycle.  For  present  purposes  a  given  cell  may  be  in  one  of 
two  states:  functionally  active,  i.e.  or  undergoing  mitosis,  i.e.,  mitotic;  other 
detailed  aspects  of  the  cycle  are  temporarily  ignored.  During  the  functionally 


active  state  occupancy  time  a  cell  performs  its  intended  function,  i.e.,  that  of 

removal  of  nutrients  and  wastes  from  contiguous  blood  flow.  Such  a  cell  has  a 

certain  natural  life  time,  the  duration  of  which  is  influenced  both  by  the 

programmed  cell  death  phenomenon  apoptosis,  but  also  by  the  occurrence  of  an 

insult,  possibly  directly  from  an  environmental  agent  such  as  a  toxic  chemical, 

but  also  from  a  metabolite;  the  latter  cause  of  death  is  necrosis.  When  a  cell  dies  a 

chemical  message  is  sent  to  a  neighboring  cell  via  the  intercellular  media 

commanding  that  neighbor  to  divide,  or  perform  mitosis,  i.e.,  go  into  S  and  Gi- 

During  the  period  of  a  cell's  mitosis  the  original  collection  of  active  cells  is 

effectively  deprived  of  two  members:  the  dead  one,  and  the  neighbor  undergoing 

mitosis.  Note  that  a  dying  cell  need  not  always  signal  the  same  neighbor  to 

become  mitotic,  but  restriction  to  that  situation  is  not  exceptionally  special  and 

does  lead  to  quite  a  simple  initial  model. 

Model  D.l.  A  Deterministic  Model  Involving  Monogamous  Cell  Pairs. 

Consider  a  monogamous  cell  pair  that  alternatively  is  in  mitosis  and  is 

functionally  active.  Suppose  {X,-,  i  =  1,  2, ... }  are  independently  and  identically 

distributed  positive  random  variables  that  represent  the  time  spent  in  mitosis  for 

the  pair,  and  [Lu,  Lu,  i  -  1,  2,  ... }  be  the  lifetimes  of  the  two  paired  cells  when 

they  are  in  the  functionally  active  state.  We  assume  the  pair's  lifetimes  are 

probabilistic,  being  mutually  independently  and  identically  distributed,  and  are 

re-sampled  after  each  mitosis  event  period  terminates.  The  pair  dies,  in  effect, 

after  lifetime  L,-  =  min(L2i,  Lu),  after  which  mitosis  begins,  eventually  terminates, 

and  a  new  pair  lifetime  begins;  so  the  process  continues.  If  Io(t)  is  the  indicator 

function  that  specifies  the  state  of  the  pair  at  t,  so 

fl  if  both  cells  of  pair  are  in  totally  differentiated  state  at  time/; 
;     [0  otherwise 


then 

pD(t)  =  P{ID(t)  =  1  I  Both  pair  members  have  just  entered  the  totally 

differentiated  state  at  t  =  0}  (3.2) 

satisfies  a  renewal  integral  equation,  see  Feller  (1966): 

PD(t)  =  l-FL{t)  +  jlPD(t-x)G(dx)f  (3.3) 


where 


and 


Ii(t)  =  P{Ll{  <  t,L2i  <t}  =  P{L  =  max(Llf,L2l-)  <  t]  (3.4) 


G(x)  =  FL(x)*Fx(x)  =  jXQFL(x-y)Fx(dx).  (3.5) 

The  argument  for  (3.3):  starting  with  both  pair  members  freshly  back  from 
mitosis  the  pair  is  either,  (i),  still  in  the  functionally  active  state  at  time  t,  an  event 
of  probability  1-Fj/f)  -the  first  right-side  term  of  (3.3)  --  or,  (ii),  the  pair  has 
cycled,  first  through  its  functionally  active  state  and  then  through  mitosis,  with 
the  pair  starting  a  new  life  at  time  x;  this  is  represented  by  the  integral  term  of 
(3.3). 
Asymptotics 

Renewal  theory  results,  cf.  Feller  (1966)  or  Asmussen  (1987),  tell  us  that,  as 
time  gets  long,  a  limit  is  approached,  expressed  as 

where  E[L]  is  the  mean  time  of  pair  residence  in  the  totally  differentiated  state, 
and  E[X]  is  the  mean  time  of  cell  pair's  absence  for  mitosis.  Thus  the  productive 
fraction  of  the  time,  pp,  depends  only  upon  the  means  of  the  state  residence 
times  and  not  upon  details  of  their  distributions.  Of  course  different  cell  types, 
and  those  with  different  locations  in  an  organ,  might  well  have  different  means, 
and  hence  different  prj-values.  Moreover,  both  means  could  be  expected  to 
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respond  to  a  toxic  chemical,  so  we  would  write  £[L;  T]  and  E[X;  T)  where  here  T 
stands  for  level  of  toxin,  past  and  present.  At  the  moment  no  attempt  is  made  to 
represent  the  dependence  of  the  above  means  on  toxin  level  in  a  mechanistic 
manner,  but  that  is  an  ultimate  objective. 
Groups  of  Cells;  (Sub)  Organs 

Suppose  there  are  C  pairs  of  cells  in  an  organ,  with  pairs  behaving  as 
discussed,  and  independently  cycling  between  states.  Then  standard  limit  theory 
suggests  that  if  D(t)  represents  the  number  of  pairs  totally  differentiated  cells 
present  in  the  organ  at  time  t, 

D[  >~    \tr\      (i  \  (      } 

is  approximately  normal /Gaussian  distributed  for  large  C.  Divide  numerator 
and  denominator  of  the  r  h  s  by  C  to  conclude  that  as  C  becomes  large  the 
concentration  of  functional  cells  D(t)/C  approaches  a  constant  (in  this  case)  i.e. 
PD- 

Model  D.2.  A  Differential  Equation  Representation  of  Cell-Toxin  Interaction. 

Let  D(t)  represent  the  random  number  of  cell  pairs  that  are  in  total 

differentiated  state  at  time  t,  and  M(t)  be  the  number  of  (potential)  pairs 

undergoing  mitosis;  D(t)  +  M(t)  =  C,  a  constant  representing  the  number  of  cell 

pairs  in  the  organ.  Suppose  the  process  {D(t),  t  >  0}  is  Markov  with  generator 

defined  by 

P{D(t  +  dt)  =  D(t)  +  l|D(f )}  =  n{T)(C  -  D{t))dt  +  o(dt) 

P{D(t  +  dt)  =  D(i)  -  1|D(0}  =  A  (T)D(t)dt  +  o{dt)  (38) 


Then  if  T  is  a  constant  level  of  toxicity  in  the  system  one  can  take  expectations 
throughout  to  obtain  a  linear  differential  equation  for  D(t)  =  E[D(tj\  the  mean  of 

D(t): 


dD 
the  solution  of  which  is 


(  =»(T)-C-(ii(T)  +  l(T))D(t);  (3.9) 

at 


D(tUD(0V{tl{T)+m)t  +    ,Cf(?,   Nfl-rfr(r)+A(r»'l  (3-10) 

v  '        K  '  /i(T)  +  A(T)l  J 

Hence  the  long-run  expected  number  of  totally  differentiated  pairs  is 

limD(f)=      C^\\   N  (3.1D 

which  bears  close  resemblance  to  Model  1  in  the  long  run,  expression  (3.6),  with 
H(T)-l  =  E[X]  and  ACT)"1  =  E[L].  For  the  above  model  it  is  again  easily  seen  that  (a 
scaled  version  of)  D(t)  is  approximately  Normal  for  large  C. 
Cell-Toxin  Interaction 

The  above  motivates  the  following  deterministic  differential  equation 
formulation.  We  allow  (3.9)  to  become 

^  =  /i(f(0)C-(/i(f(f))  +  A(f(0))D(0  (3.12) 

^(Q     ,f)     P(Qf(0 

dt  u  l+KT(r) 
The  differential  equation  for  T(t),  roughly  interpretable  as  the  mean 
concentration  of  toxin  present  at  f,  expresses  the  fact  that  the  increase  in  toxin 
present  at  t  equals  the  rate  of  input  of  toxin  minus  the  rate  of  output;  the  latter 
rate  is  modeled  as  a  saturated  (Michaelis-Menten  type)  service  rate, 
vT(t)/(l+  KT(t)),  multiplied  by  the  concentration  of  totally  differentiated  cells 

in  the  fluid.  In  a  later  section  we  study  a  complete  probabilistic  version  of  (3.12). 
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It  is  perhaps  natural,  but  is  conjectural,  to  attribute  functional  dependence  of 

the  cell-change  rates  in  a  way  that  tends  to  reduce  the  rate  of  increase  of 

functional  cells,  the  "toxin  fighters",  with  toxin  concentration.  For  example 

(only),  let 

H(T)  =  »0e-ST;    A(r)  =  A0^T  (3.13) 

where  £  and  rj  are  positive.  This  choice  tends  to  kill  more  cells  as  T  increases  but 

to  increase  the  number  undergoing  mitosis  since  the  duration  of  the  mitotic 

period  is  increased.  Another  possibility  is  that  the  condition  of  cells  emerging 

from  mitosis  is  changed  in  the  presence  of  a  toxin:  some  such  cells  are  deficient, 

perhaps  leading  to  birth  defects,  or  others  become  initiated  for  cancer.  We  do  not 

now  model  these  many  options.  Supposing  r(t)  =  r,  a  constant,  the  long-run  toxin 

concentration  is  the  solution  of 


D-S^) 


and 


v 


DT 


or 


r  = 


KA(T) 

vCn(T)T 
1+kT  "(l+Kf)[/i(f)  +  A(f)] 

vCn0T 


(3.14) 


(3.15) 


A  graph  shows  that  a  finite  steady-state  long-run  solution  for  toxin  level,  namely 
f,  will  exist  under  certain  conditions.  Note  that  in  Fig.  1  two  solutions  to  (3.15) 
can  exist. 
If 


( 


x  <  max 

7 


> 


vCjd0T 


(i+kFJuo  +  V^1 


j; 


Right-hand  side 


Left-hand  side 


►    T 


T 


Figure  1 
It  seems  clear  that  the  smaller  of  the  two  possible  solutions  is  biologically 
plausible,  since  it  is  evidently  an  increasing  function  of  t,  whereas  the  other 
solution  decreases.  If  the  above  inequality  is  not  satisfied  the  organ  is  soon 
overcome  and  toxin  concentration  increases  indefinitely. 

Model  D.3.  A  Cell-Age-Specific  Differential  Equation  Model  of  Cell-Toxin 
Interaction:  The  Effects  of  Apoptosis  and  Necrosis. 

In  order  to  model  the  effects  of  programmed  cell  death  or  apoptosis,  and  death 
from  insult,  called  necrosis,  it  is  necessary  to  represent  the  cell-aging 
phenomenon.  We  ch  Jose  to  do  this  in  the  context  of  Markov  chains  or  rate 
equations,  simply  extending  (3.12),  by  the  consideration  of  two  totally 
differentiated  cell  types:  new  or  young,  namely  those  that  have  comparatively 
recently  been  "born",  i.e.,  emerged  from  mitosis,  and  old,  those  that  have 
transitioned  from  the  new  state  ("reached  maturity")  and  are  now  eligible  to  die. 

The  following  rate  equations  reflect  the  effects  of  a  toxin  in  the  above  process. 
Let  Dn(t)  be  the  number  of  new  or  young  cell  pairs  (recently  emerged  from 
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mitosis),  Do(t)  the  number  of  old  cell  pairs,  and  T(t)  the  amount  of  toxin  in  contact 
with  them.  Then  we  write 

^l  =  -A„(e^(O-l)D„(()-0Dn(f) 

+^-mO(C-D„(()-D0(<))  (316) 

^l  =  -Xoe^n')Do(t)  +  0D„(t)  (3.17) 

<^(<)_Tm     v„f(f)Dr(f)     t)0f(/)D0(() 
d»         ()-   l+^f(0   "   1+^(0 

Consider  first  (3.16).  The  first  term  on  the  right-hand  side  represents  a  death 
rate  of  new  cell  pairs  that  increases  with  toxin  concentration  and  is  zero  with 
none.  It  represents  the  effect  of  toxin-induced  necrosis  upon  the  new  cell  pairs. 
The  second  term  is  an  aging  term:  it  represents  the  rate  at  which  new  cells  enter 
the  old-cell  population.  The  third  term  represents  the  mitosis  rate;  mitosis 
resembles  birth  in  that  new  cell  pairs  are  the  product.  It  has  been  assumed  that 
mitosis  rate  always  decreases  with  increased  toxin  concentration,  but  such  need 
not  be  the  case. 

Next,  consider  (3.17).  The  first  term  is  a  death  rate  for  old  cell  pairs  that 
increases  with  toxin  concentration;  it  represents  both  necrosis  and  apoptosis.  The 
second  term  is  the  rate  of  maturation  of  new  cells  into  old.  Doubtless  the 
transition  rate,  0,  could  depend  upon  toxicity  in  the  organ,  but  this  possible 
dependence  is  not  modeled  here.  Biological  insight  and  information  can  perhaps 
suggest  an  appropriate  and  interesting  functional  dependence  for  <p. 

Finally,  (3.18)  models  the  rate  of  toxin  concentration  increase.  The  first  term  is 
the  rate  at  which  the  basic  toxin  reaches  the  vicinity  of  the  cells.  The  second  and 
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third  terms  both  represent  the  rates  at  which  toxin  concentration  is  reduced  by 
action  of  new  cells  (second  term)  and  old  cells  (third  term).  Note  that  bom  cell 
typ  ■  reduce  toxin  concentration  roughly  in  proportion  to  that  concentration,  but 
their  effects  saturate  a  la  Michaelis-Menton. 

The  equations  supplied  are  illustrative  of  what  may  happen  for  certain  toxins 
and  cell  types;  they  are  intended  to  provoke  discussion  and  stimulate 
experimental  and  observational  checking  and  revision. 

A  steady-state  solution  to  equations  (3.16)  -  (3.18)  for  a  constant  toxin  input 
t(/)  =  t  would  satisfy 

0  =  -Xn(e^T  -l)Dn-<pDn+iie--T(C-Dn-D0)  (3.16a) 

0  =  -  V^Dq  +  0Dn  (3.17a) 

0  =  T-i^L-"0TP0  (3.18a) 

1  +  KnT      1  +  KqT 

Simplification  of  (3.16a)  -  (3.17a)  results  in 

DQ=4-e-^TD„  (3.19) 


A 


0 


Dn=^0e   ^Cx 


An(e7?"T-l)  +  0  +  ^"^T 


(3.20) 


Finally,  expressions  (3.19)  -  (3.20)  can  be  substituted  into  equation  (3.18a)  to 
obtain  an  equation  for  T.  This  latter  equation  can  have  0, 1  or  2  solutions.  If  it  has 
no  solutions,  then  Dn  =  Dq  =  0  and  the  organ  is  dead.  If  the  equation  has  two 
solutions,  then  the  smaller  of  the  two  possible  solutions  Ts  is  biologically 
plausible  since  the  smaller  solution  is  an  increasing  function  of  t,  whereas  the 
larger  solution  is  a  decreasing  function. 
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If  t  =  0,  then  T  =  0  and 


Dn= p  (3.21) 


Dn  = 


(3.22) 


Mui). 


4.  Numerical  Illustrations. 

We  have  submitted  the  above  hypothetical  equations  to  several  dosage 
regimes. 
Illustration  1. 

Suppose  toxin  input  is  0  until  time  t  =  5,  after  which  a  constant  exposure  r  is 
given.  Figure  1  displays  the  function 

f(T)=  VnTDn  i  V°TD° 

}K   '     \+KnT     l+K0T 

where  Dn  and  Do  are  obtained  from  equation  (3.19)  and  (3.20);  the  parameter 
values  used  are  displayed  on  the  figure.  The  (possible)  solution(s)  to  equation 
(3.18a)  can  be  found  drawing  a  horizontal  line  at  the  value  of  r  of  interest.  The 
limiting  amount  of  toxin  in  the  organ  subject  to  a  constant  input  of  toxin  t  is  the 
smaller  of  the  two  values  of  T  corresponding  to  the  intersection  points  of /with 
the  horizontal  line.  Note  that  there  is  no  solution  for  x  >  3.53,  the  maximum  value 
of  the  function.  This  means  that  toxin  has  killed  all  old  and  new  cells;  long  before 
that  occurs  the  organism  has  died. 

Figures  2-5  display  the  results  of  a  numerical  solution  of  equation  (3.16)- 
(3.18)  for  values  of  t  =  0, 1,  2,  3,  3.3,  and  3.53.  The  other  parameter  values  used 
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are  displayed  on  the  Figures.  The  initial  values  of  the  numbers  of  young  and  old 
cell  pairs  and  the  amount  of  toxin  present  are  the  limiting  values  when  x-  0. 

Figure  2  shows  the  effects  of  various  levels  of  toxin  on  the  young  cell  pair 
population.  At  the  beginning  of  exposure  the  young  cell  population  decreases. 
However,  if  the  level  of  toxin  is  not  too  large,  after  this  initial  decrease  the  young 
cell  population  increases  to  a  limiting  value  higher  than  that  occurring  when 
r  =  0.  If  the  level  of  exposure  is  too  high,  the  young  cell  population  decreases  to 
zero. 

Figure  3  shows  the  effects  of  various  levels  of  toxin  on  the  old  cell  pair 
population.  If  the  level  of  exposure  is  too  high,  the  old  cell  population  decreases 
to  0.  For  low  or  moderate  exposure  levels  the  old  cell  population  initially 
decreases  but  then  recovers  somewhat  to  a  limiting  value  which  is  below  the 
limiting  value  if  there  were  no  toxin. 

Figure  4  displays  the  total  number  of  totally  differentiated  cell  pairs.  As  the 
level  of  exposure  increases  the  total  number  of  cell  pairs  decreases.  The  total 
number  of  totally  differentiated  cell  pairs  initially  decreases,  then  recovers 
somewhat  before  decreasing  a  little  again  to  a  limiting  value. 

Figure  5  displays  the  level  of  toxin  in  the  organ. 

It  seems  possible  that  our  model  represents  an  effect  called  hormesis  wherein 
a  smallish  amount  of  deleterious  agent  induces  a  biological  entity  to  over- 
compensate  for  a  (low-level)  insult,  but  eventually  to  succumb  to  a  larger  dose.  If 
young  cells  are  actually  more  efficient  and  productive  than  old  then  the 
proportional  buildup  of  the  new-cell  population  for  low  toxin  levels  has  just  such 
an  effect  in  our  example. 
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5.   Stochastic  Versions  of  Models 

The  purpose  of  this  section  is  to  describe  a  fully  stochastic  version  of  the  cell- 
toxin  interaction.  Note  that  once  stochastic  elements  are  permitted  to  enter  a 
great  many  optional  behaviors  can  be  modeled.  It  will  be  seen  that  expressions 
(3.12),  (3.16),  (3.17),  and  (3.18)  are  not  exact  replicas  of  the  correct  mean  or 
expected-value  equations. 
Model  S.l.  Markovian  Model  of  Cell-Toxin  Interaction. 

Suppose  D(t)  represents  the  number  of  totally  differentiated  cell  pairs  in  an 
organ  with  a  maximum  of  C  pairs,  and  T(t)  is  the  toxin  concentration  at  time  t. 
Both  D(0  and  T(t)  are  now  considered  to  vary  randomly. 

Let  (D(t),  T(t),  r  >  0)  be  a  bivariate  birth  and  death  Markov  process  in 
continuous  time  and  with  state  space  Rx  R,  i.e.,  all  pairs  of  integers  such  that 
0  <  D(0  <  C,  0  <  T(t).  Then  put  for  the  transition  probabilities  from  t ->  t  +  A 

P{{D(t  +  A),T(*  +  A))  =  (D(0  +  1,T(0)|D(0,T(0}  =  »{T(t)){C  -  D(t))A  +  o(A) 

P{(D(t  +  A),T(*  +  A))  =  (D(0  -  l,r(f))|D(0,T(0}  =  A(T(f))D(0A  +  o(A) 

P{{D(t  +  A),r(f  +  A))  =  (D(0,T(0  +  1)|D(0,T(0}  =  r(t)A  +  o(A) 

P{(D(f  4-  A),  T(f  +  A))  =  (D(t),  1(0  -  l)|D(f ),  r(f)}  =  ^^J^  ^  +  o(A) 

P{(D(f  +  A),  T(*  +  A))  =  (D{t),T(t))\D(t),T(t)}  =  1  -  p{D(t),T(t))A  +  o(A),        (5  1} 

where  p(D(0,T(0)  =  /i(T(0)(C -  D(f ))  +  A (T(f ))D(f )  +  <f)  +  ^(0^(0  /  (l  +  *T(0); 
1-pA  is  simply  the  probability  of  no  change  in  the  time  period  (t,  t  +  A). 

Examination  of  (5.1)  reveals  that,  under  very  special  conditions,  expected- 
value  and  higher-moment  equations  can  be  written  down,  and  can  be  solved 
numerically.  Non-linearities  in  general  impede  such  a  step,  but  approximation 
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and    perturbation    methodologies    can    be    invoked    to    produce    useful 
approximations. 
The  Mean  Equations 
Invoke  (5.1)  to  obtain 

E[D(t  +  A)|D(0,T(0]  =  (D(0  +  l)/i(T(0)(C  -  D(t))A 

+(D(0-1)A(T(0)D(0A 

4  D(0[l  -/i(T(0)(C  -  D(0)A  -  A(T(f))D«A]  +  o(A) 

=  fi{T(t)){C  -  D(0)A  -  A  (T(O)D(OA  +  D(f)  +  c(A).       (5  2) 

Hence,  upon  removing  the  condition  and  passing  to  the  limit 


dt 


E[D(0]  =  E{v(T(t))]C  -  £[D(0/i(T(0)]  -  E[A(T(0)D(f)] 


(5.3) 


Next, 


:[r(f  +  A)|D(0,r(0]  =  (T(0+i)r(OA  +  (r(0-i)t,^MA 


U[      t  l+*T(f)    Jj 


Removal  of  the  condition  leads  to 


(5.4) 

These  resemble  the  previously-presented  deterministic  model  (3.12),  but  only 
roughly,  since  expectations  of  non-linear  functions  such  as  appear  above  cannot 
be  directly  evaluated.  Special  methods  can  be  used  that  give  satisfactory 
approximations. 
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Clearly  the  equations  (5.3)  and  (5.4)  cannot  be  solved  without  the  introduction  of 
approximate  equations  for  the  non-linear  expectations  on  the  r  h  s.  One  approach 
for  doing  this  is  presented  in  Appendix  A.  A  second  approach  is  presented  in  the 
next  section. 
Asymptotics  for  C  Large. 

Since  the  number  of  cells  in  an  organ,  C,  is  large,  being  of  order  107,  it  is 
natural  to  consider  explicit  asymptotic  approximations  based  on  that  fact.  In 
what  follows  we  proceed  formally.  For  precedent  here  see  McNeil  and  Schach 
(1973)  but  also  more  recent  work. 

Define  the  joint  moment  generating  function  (assumed  to  exist),  to  be 


v(ed,er,t)  =  Ele9dD{t)+mt) 

Use  the  Markov  property  to  derive  forward  equations  in  transform  space: 

e^(D{o+i)+^T(0/i(T(0)(c_D(f))A+e^(D(0-i)+etr(0A(T(0)D(OA 

eMthe,(T(th\)  ,  u     frD(0+gt(r(Q-i)  vD(t)T(t)A 

K)  1+kT(0 

+e  W>+ W)[!  _  p(D(t),T(t))A]  +  o(A), 


where 


p(D(t),  T(t))  =  M(T(f))(C  -  D(0)  +  A  (T(O)D(O  +  z{t)  +  ^™- 


(5.5) 


17 


Thus 

'e6^D^+d^e^%{T(t))(C-D(t))A 


=  E 

+E 

+E 
+E 


gWhW%(r(t))D(t)A 

E^D(0+W,)p(DW/T(0)J 


,edD(t)+otT(t) 


Consequently,  after  re-arrangement  and  division  by  A,  and  letting  A  — » 0, 

mSidH^i  -yH^('^(r(i))(c-D(f))| 

,ewD(fHe>nn  p(0r(0 

1+kT(0 


■u(f)(e"e'  -l)E 


(5.7) 


These  equations  are  highly  non-linear  and  intractable;  consequently  we  introduce 

an  appropriate  asymptotic  analysis  that  assumes  C,  the  number  of  cell  pairs,  to  be 

large. 

Scaling 

Define 


x(0 


_  D(Q-C«W   y(0  =  T(t)-Cp(t) 


Vc      ,lv>~      Vc 

and  the  joint  mgf  of  the  scaled  "noise"  variables  X(i)  and  Y(t): 


(5.8) 
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(d)    r 
<p{9d,0t;t)  =  E 


,edx(t)+etY(t) 


=  Er  edD(t)/Jc+BtT(t)/Jc~l -(ed-Jca(t)+ety[cp(tj) 


(5.9) 


Consequently, 


^{ed/^tej^^  =  cp(edA^[6Mi)+m)y  (s.io) 

It  is  necessary  to  find  equations  for  a(t)  and  p(t),  and  the  joint  mgf  (p{6d,  6f,  t). 

Proceed  as  follows: 

use  (5.7);  re-define  the  transform  variables  as  below  6d  -»  6d  /  Vc,  6t  — »  6t  I VC ; 

from  (5.10)  we  get 

d*(ed/  Vc,  et  I  Vc;f)  JJL^M^m)  +  ^^(t) + ^'(0)^(0-a(0+fl|/,(0) 


=  E 


d/ 


^^/VcJtVcxio+CaCoy^/VcXVcYCO+c/Jio) 


x{(> /V^  -  l)(/i(Q3(0  +  VC Y(0)(C -  Ca(0  -  VCX(* 


X-^/Vc  _ ! VA (C^)  +  VCY(0)(Ca(0  -  VCX(* 
^/Vc.^ 

wf/e-*,/Vc     N(ca(0^Vcx(0)(cw)+Vcy(0)^ 

UV  J  1+K:(Q3(0  +  VCY(0) 

Next  scale  the  transition  rates 

A*(c«0+VcY(0)=/(«0+y(0/Vc) 
a(c«o+Vcy(0)  =  a*(/3(0+y(0/Vc) 


(5.11) 


(5.12a) 
(5.12b) 
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v(t)  =  v(t)/C  <t)  =  Cz\t)  (5.12c) 

K(cp(t)  +  VCY(O)  =  k  (p(t)  +  Y(t)  I VC)  (5.12d) 

in  order  to  obtain  non-degenerate  results  involving  all  aspects  of  the  process  as  C 
becomes  large.  Starred  functions  or  parameters  are  0(1). 

Expand  state-dependent  parameters  by  Taylor  series  in  powers  of  1  /  VC  after 
cancelling  ^(^a0+ W))  from  both  sides;  we  put  n*Q(p(t))  for  the  first  term  of 

the  expansion  of  //(/J(f)  +  Y(f)/Vc);  notation  for  A*  and  k*  is  analogous.  We 
obtain 


,B,X{t),e,Y[t) 


\ 


+ 


ed    .  ed  . 


C     2C 


-^+^+...Y/iSWf))+rt(W))-^-)c-Ca(f)-Vcx(f)) 

...Y^(«0)+AiWO)-^-)c«(0+-^x(0) 


+  — *-+... 


VC     2C 


c/(0 


y 


+ 


et     ef 

VC     2C 


u*(0/c 


(co(0 + Vcx(0)(c^(f ) + Vcy(O)' 


i+^08(O)+»rlO(O)y(O/Vc 

Re- write  the  last  line  above  in  expanded  form  to  get 


(5.13) 


l  +  *o(/»(0) 


.-(«iW0)/(i+«5W0))ffi 


\ 


:  + 


(5.14) 
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&+Jc(eia'(t)+etp'(t))<p= 


u6d    ej     ^ 

'   +-B-  + 


VC      2C 


c 


J 


1    dcp 


/ij(0(i  -  «(*))* + nlW  ~  «(0)^  J^  -  /i5  (0 


V 


1     do        1 


'C<?0 


VC     2C 


c 


Ao(()a(0,  +  A-l(()«(0^^  +  Ao«)^^+...) 


+  — L-+.. 


VC     2C 


C<px\t) 


J 


JC     2C       , 


C 


»,  ,a(tWt) 
1+jc0(0 


(( 


1      <fy 


ICddr. 


+  . 


+;w4Lii.;(()«)^,„ 


(5.15) 


We  express  the  dependence  of  <p  on  C  as  follows: 


(5.16) 
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Now  substitute  into  (5.15)  to  obtain 


0{w1+^a'{,)+e<P'%*{lc]-- 


;=0 

f 


ed    ,  ed   , 


VC     2C 


c 


y 


,o(0(l-«(0)I;;{^)+,;(0(l-a(0)^I^— 


^cj=oded  Vc; 


V 


C      2C 


y 


1    $d<Pj      1 


/To  ;(Vc);  ^;4^(VC) 


1    $d<Pj      1 


vc/=o^(Vc); 


( et     e? 


VC     2C 


> 


/=o    (vc) 


+ 


VC     2C 


c 


ui+Kb(0/tS  '(Vc)' 


*o(0  icpbded{jcy 


v(t)a(t) 


l+KQ(t) 


1      1(0*1(0 

i+k5(0 


M 


1     y  gg;        1 

^r;ts^(vcy+"j> 


where  the  omitted  terms  ("...")  are  0(1 /C),  and  Aq(0  =  ^0  WO)  etc. 


(5.17) 
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Identify  terms  of  order  Vc : 
lhs:(eda'(t)  +  etp'(t))<pQ 

r  h  s:  ed[(So(p(t))(l  -  a(t))) -  a(t)^0(p{t))}l>o  +  *. 


T*(t)- 


V(t)a(t)p(t) 


(5.18a) 


(p0;  (5.18b) 


i  +  ^o(«0) 

since  these  must  cancel,  the  functions  exit)  and  Pit)  satisfy  the  differential 
equations 


«'(0  =  ^W0)(l-a(0)-Ao(«0)a(0 

a(t)0(t) 


P'(t)=t*(t)-v*(t) 


i+*o(«0) 


(5.19) 


the  latter  must  in  general  be  solved  numerically,  but  steady-state  information  can 
be  derived  more  easily.  The  system  of  equations  (5.19)  is  a  scaled  version  of 
equation  (3.12). 

Next,  identify  terms  of  order  1  (coefficient  of  (l/Vc)  .  The  expression 

obtained  is  of  the  form 


dt  2     dW      2     nj 


V°J 


+B**m^+Bdtm^ 


+Btd(t)0d^+*ttm 


(5.20) 
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Here  the  above  coefficients  are 

4*(0    =  (/io(0(i-«(0)+4(0«(0) 
{  i+«b(0   J 

%(0    =  /*i(0(i-«(0Mi(0«(0 


Btt(t)      =  ~ 


v(t)a(t) 


1  g(0a(0 


(5.21a) 
(5.21b) 

(5.21c) 
(5.21d) 

(5.21e) 

(5.210 


i+ir0(0 

The  resulting  partial  differential  equation  for  (f>Q  is  recognizable  as  characterizing 
the  joint  mgf  of  an  Omstein-Uhlenbeck  (Gaussian)  stochastic  process. 
Moment  Equations. 

Differentiation  of  the  pde  with  respect  to  $d  and  0/  allows  recovery  of 
differential  equations  for  the  covariance  function  of  (X(f),  Y(t)),  the  scaled 
stochastic  terms  describing  deviations  from  the  mean  (to  order  C)  term 

E[D(0]  =  Ca(t),  E[T(t)]  =  Cp(t)  (5.22) 

Differentiation  of  <Po{6d,et;t)  at  6d  =  6t=0  shows  first  that  if  X(0)  =  Y(0),  then 
E[X(t)]  =  E[Y(t)]  =  0  to  the  order  suggested.  Second  derivatives  at  6  =  0  then 
deliver  these  differential  equations  for  the  variances  and  covariances, 
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Var[X(t)]  =  E[x2(0]  =  mdd(t),Var[Y(t)}  h  e[y2(0]  =  mtt(t)  (5.23) 

Cov[X(t),  Y(t)]  =  E[X(t)  ■  Y(0]  =  ma  (0,  (5.24) 

jtmdd(t)  =  Ad(t)  +  2Bdd(t)mdd(t)  +  2Btdmdt(t)  (5.25) 

^^=(B^(0  +  B«(0W(0  +  Bj/(0^rf(0  +  B^«(0         (5.26) 
^  m«  (0  =  ^  (0  +  2BW  (0««  (0  +  2B^mrff  (0-  (5.27) 

These  equations  can  be  solved  numerically.  Owing  to  the  appearance  of  the 
Ornstein-Uhlenbeck,  Feller  (1966),  we  can  argue  that  for  large  C,D(t)  is 
approximately  normal  with  mean  Ca(t)  and  variance  C  mdd{t)',  T(t)  is 
asymptotically  normal  with  mean  C  p(t)  and  variance  C  m^(f);  the  two  quantities 
are  correlated  with  correlation  coefficient  pdt (t)  =  mdt{t)  I  yjmdd(t)mtt(t). 

Steady-State  Solutions:  Dose-Response  for  Steady  Exposure 

Suppose  r*(0  =  x ,  a  constant,  and  that  exposure  has  proceeded  for  some  time 
so  that  a  steady-state  condition  has  been  reached;  this  is  modeled  by  letting 
0^(0  =  p(t)  =  0  in  (5.19).  We  have,  suppressing  t  and  invoking  the  tentative 
working  parametric  forms  (3.13), 

/i0r^(l  -  a)  =  V^a  (5.28a) 

T*  =  t/__«L_.  (5.28b) 

l  +  K0{p) 

Notice  that  (5.28a)  can  be  solved  for  /?,  the  steady-state  toxin  concentration,  in 

terms  of  a,  the  steady-state  fraction  of  totally  differentiated  cells: 


Vo 


(l-aW 


^0 


(5.29) 


V    or    J) 

This  can  now  be  inserted  into  (5.28b)  to  obtain  an  implicit  dose-response 
relationship: 
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T    =  V 


a 


In 


Wi-<* 


(5.30) 


i+*o(W({* +»f)   UoV  «  ;/ 

the  above,  (5.30),  can  be  solved  for  a[x*]  so  as  to  depict  the  fraction  of  totally 

differentiated  cells'  dependence  on  input  toxin;  the  latter  must  here  be  a  steady 

flow. 

Slope  of  Dose-Response  Curve  for  Small  Dose. 

Of  interest  to  risk  analysis  is  the  behavior  of  the  dose  response  curve  for  small 
values  of  (toxic)  dose.  We  approximate  this  by  finding  an  expression  for 
d 


da 


dx 


=   — r(l-tt  T    j    =   '  tne  rate  °f  increase  of  the  fraction  of  non- 


functioning  totally  differentiated  cells  when  toxin  input  is  very  small. 

First  notice  that  if  i  -  0  then  /i0(l  -  a)  I  A0a  =  1 .  So  differentiation  of  (5.30)  at 
T*  =  0  provides 


da 


dx 


T  =0 


ft*  *\ 

i  +ri 

» 

V 


{l+4(pp-a[0]) 


Vo  +  h 


(5.31) 


Another  differentiation  will  give  information  concerning  the  curvatures  of  the 
response,  or  propensity  to  exhibit  "hockey-stick"  or  knee-like  or  threshold 
behavior  at  small  dose  levels.  A  knowledge  of  such  behavior  is  of  interest  to 
toxicologists  and  regulators  who  are  concerned  with  risks  in  the  workplace.  It 
will  be  kept  in  mind  that  trie  behavior  elucidated  depends  to  an  unknown  degree 
upon  the  suitability  of  the  models. 
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Model  S.2.  Markovian  Model  of  Cell- Age-Specific  Cell-Toxin  Interaction. 

Suppose  Dn(t)  represents  the  number  of  young  (or  new)  cell  pairs  (recently 
emerged  from  mitosis),  D0(t)  the  number  of  old  cell  pairs,  and  T(0  is  the  toxin 
concentration  at  time  t.  All  the  variables  Dn(t),  D0(t)  and  T(t)  are  now  considered 
to  vary  randomly. 

Let  (D0(t),  Dn(t),  T(t);  f  >  0)  be  a  trivariate  birth  and  death  process  in 

continuous  time  and  with  state  space  Rx  Rx  R,  i.e.,  all  triples  of  integers  such 

that  0  <  D0(t)  +  D„(f)  <  C,  0  <  T(t).  Then,  put  for  the  transition  probabilities  from  t 

-*t  +  A. 

PftD^f  +  Ap^f  +  AjT^  +  A^fD^O  +  l^D^O-lTWJiDo^D^OTW} 

=  0Dn(OA  +  o(A) 

P{{D0(t  +  A),D„(f  +  A),T(f  +  A))  =  (Do(0  -  l,Dn(0,T(0)|Do(0,Dn(0,T(0} 

=  A0(T(f))D0(f)A  +  o(A) 

P{{D0(t  +  A),DB(f  +  A),T(f  +  A))  =  {D0(t),Dn(t)  +  l,T(0)|Do(0,Dn(f),T(0} 

=  /i(T(0)(C-D„W-Do(0)A  +  o(A) 

P{(Do(f  +  A)/Dn(f  +  A)/T(f  +  A))  =  (Do(0,D„(0-l,T(0)|Do(0,Dll(f)/T(0} 

=  Afi(T(f))Dn(f)A  +  o(A) 

P{(D0(f  +  A),DB(f  +  A),T(t  +  A))  =  (D0(t),Dn(t)J(t)  + l)N0^n(0^«} 

=  r(0A  +  o(A) 

P{(Do(f  +  A)/Dn(f  +  A),T(f  +  A))  =  (Do(0,Dn(0,T(0-l)|Do(0,D„(0/T(0} 

1  +  k-„T(0  l  +  r0T(f) 

P{{D0(t  +  A),DR(f  +  A),T(*  +  A))  =  {D0(t),DMnt)Po(t).Dn(*)>T(t)} 
=  l-p(Do(0,Dn(f),T(0)A  +  o(A) 
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where 


KDo(0,Dn(0/T(0)  =  ^Dn(0  +  Ao(T(0)Do(0+^(TW)(C-Dn(0-Do(0) 
+An(T(0)D„(0  +  T(0 

+  ^no  D(0+_j^LD(0 

l+rMT(0    nW    1+jcoT(0    oW 
Asymptotics  for  C  Large. 

In  this  subsection  we  invoke  the  techniques  applied  for  the  Markovian  model 
of  Cell-Toxin  interaction  to  obtain  a  limiting  trivariate  normal  distribution  for  the 
number  of  young  and  old  cell  pairs  and  amount  of  toxin  in  the  organ  as  the 
number  of  cells  in  the  organ  C  — > «». 

Define  the  joint  moment  generating  function,  (assumed  to  exist),  to  be 

y{e0,en,et-,i)  =  E[exP{e0D0(t)+enDn{t)+etT(t)}].  (5.32) 

Use  the  Markov  property  to  derive  the  forward  equations  in  transform  space: 
E[exp{e0D0  {t  +  A)  +  enDn(t  +  A)  +  9tT(t+  A)p0  (t ),  Dn  (t),  T(tj\ 
=  exp{eo{Do(t)  +  l)  +  6n{Dn(t)-l)  +  eir(t)}0Dn(t)A 
+exp{0o(Do(O-l)  +  0nDn(f)  +  ^T(r)}Ao(T(O)Do(OA 
+exp{0oDo(O  +  0n(Dn(O  +  l)  +  ^T(O}M(T(O)(C-Dn(O-Do(O)A 
+exp{0oDo(O  +  0„(Dn(O-l)+^T(O}An(T(O)Dn(OA 
+exp{0oDo(t)  +  enDn(t)  +  et{T(t)  +  l)}r{t)A 

+exp{60D0(t)  +  enDn(t)  +  9t(T(t)-lj\ 

+  exp{tf0D0  (0  +  6nDn  (0  +  0,  (T(f ))}[!  -  p(D0  (f),  D„  (f),  T(0)a]  +  o(  A). 


i;nr(f)Dn(Q  (  vJ(t)D0(t) 
I  l  +  KnT(t)        l  +  K0T(t)  _ 


(5.33) 
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After  rearrangement  and  division  by  A,  letting  A  ->  0,  and  putting  Z(t)  = 

exp{e0D0(t)  +  enDn(t)  +  etT(t)} 

dt 
=  (eeo-e»-l)E[Z(t)0Dn(t)] 

+{e-0° -l)E[Z(t)X0(T(t))D0(t)] 

+{ee»  - l)4z(0/i(T(f))(C -  D0(t) -  Dn{t))] 

+(e-e"-l)E[z(0^(T(0)Dn(0] 

+(e9< -l)t(t)E[Z(t)] 


+(*-*' -l)E 


Z(t) 


vnT(t)Dn(t)     v0T(t)D0{t) 


+ 


l  +  KnT(t)  l  +  K0T(t) 


Define  as  before 


(5.34) 


MrD0(f)-Ctt0(Q^       (f)rDB(0-C«„(QJ(f)  =  T(0-Ci3(0 

Vc  Vc  Vc 

and  the  joint  mgf  of  the  scaled  "noise"  variables  X0(t),  Xn(t),  Y(t): 

9{o0fen/et-,t)  =  E[exP{e0x0(t)  +  enxn(t)  +  etY(t)}} 


(5.35) 


(5.36) 
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Proceeding  as  in  the  derivation  of  (5.11) 
%exp[JC(d0a0(t)+enan(t)+etp(t)j\ 


dt 


+(fylc{60ai(t)  +  6na'n  (t)  +  etp'(t))exp{*Jc(60a0(t)  +  6nan  (t)  +  6tp(t)j\ 


=  E 


exp{(0o/Vc)(VcXo(r)  +  Cao(O)  +  (^/Vc)(Vcxn(O  +  Ca„(O) 


+(^/Vc)(Vcr(0  +  c«0)} 


+^0 / VC  _  iJAo(C/?(0  +  VCY(0)(Cao(0  +  Vcxo(0) 

+^-/^-l^(C«0  +  VCY(0)(c(l-ao(0-aB(0)->/CXo(0 
+(<,-*„  /  VC  _  ^  (CW)  +  VC  Y(0)(Ca„  (0  +  VCXn  (0) 


+(V*</^-i) 


un(c^(0+VcY(0)(c«n(0+^xn(Q) 
i  +  ^(cw)+Vcy(o) 

v0(C/?(Q  +  VCY(Q)(Cao(0  -f  VCX0(Q)" 

i+iro(c/5(0+Vcy(0) 


(5.37) 
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Next  scale  the  transition  rates  as  before 


/i(c«0+VcY(0)=MSWO)+^(«0)^ 

;  VC  (5.38a) 

k0(cp(t)  +  VCY(O)  =  A;,o(^(0)  +  *o,l  WO)^  (538b) 

VC  (5.38c) 


u=u  /C  r(f)  =  Cr  (0 


K"0  —  L  K"0  K"R  —  L  Kn 


(5.38d) 


(5.38e) 

Substituting  into  equation  (5.37)  and  identifying  terms  results  in  the  following 
equations. 

The  terms  of  order  Vc  imply  that  a0(t),  an(t),  and  Pit)  satisfy  the  differential 
equations 

a'o(t)  =  (pan(t)-^o0{p{t))ao(t)  (5.39a) 

«;(O  =  ^oWO)a-«o(O-«n(O)-0«n(O-^oWOK(O  (539b) 

^'M  =  T* W  ~  ■,      l*„,s  v*nPM<*n (0  "  ,      \Rt  ,  y0P(t)a0(t)-         (539c) 
1  +  ^(0  l+fo«0 

These  equations  are  scaled  versions  of  equations  (3.16)  -  (3.18). 
The  terms  of  order  (VC )    result  in  the  following  equations. 
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dt 


( <a 


^-A0(t)^An(t)^At(tye0en4>an 


\ 


<Po 


+e, 


+er 


R  d(P0     ,    R  ^0     ,R 


u 


^<p0 


-0* 


^0     n"lden     ""  ddt 


where 


i4„(0  =  /ij(«0)[l-«0(0-«n(0]  +  ^#o(«0)«n(0+^»(0 

B0,„  =  0 

Bo,f=-Aa,l03(0)«o(0 

B„/O=-^(/?(0) 

Bn,n=-^oWO)-4o(«O)-0 

B„,<  =  lii  (/3(f ))(!  -  a0  (0  -  «„  (0) -  «n  (04,1  WO) 


=  -^(0    r    .  -^(0 


Bt,t  = 


_ZJWl(0 


i+^(0 L    i+^WJ  i+f0W) 


voao(0 


1     ^(0 

1+jcJ5(0J 


(5.40) 

(5.41a) 
(5.41b) 

(5.41c) 

(5.41d) 
(5.41e) 

(5.41f) 
(5.41g) 
(5.41h) 

(5.41i) 

(5.41J) 
(5.41k) 
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Finally,  calculations  similar  to  those  leading  to  (5.23)  -  (5.27)  result  in 


dt 


fco2(0]  =  A> 


+2B0/0E 


X20(t) 

+2B0fnE[X0(t)Xn(t)] 

+2B0rtE[X0(t)T(t)} 
j-E[x2n{t)\  =  An       +2Bn>0E[X0(t)Xn(t)] 

+2Bn,„E[xn2(0] 
+2BnttE[Xn(t)T(t)] 


7tE 


T2(t)]  =  At       +2Btf0E[X0(t)T(t)] 
+2BttnE[Xn(t)T(t)] 

+2BtttE\T2(t) 


jtE[X0(t)T(t)]  =  [B0>0  +  Bl>t]E[X0(t)T(t)] 


(5.42a) 


(5.42b) 


(5.42c) 


+BO/nE[X„(0T(0]  +  BME[Xo(f)Xn(0]  +  BO/fE[T2(0] 
X2« 


+Zt>0E 


^E[Xn(0T(0]  =  Bn#oE[Xo(f)T(f)]  +  (Bn#„  +  Bu)E[XII(f)T(0] 

+B^oE[Xo(0Xn(0]  +  B^E[xn2(0]  +  BW/fE[72(0 


(5.42d) 


(5.42e) 
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^E[X0(f)Xn(f)]  =  [B0f0  +  BB#n)E[X0(f  )X„(0] 

+Bn>tE[X0(t)T(t)] 

+B0>tE[Xn(t)T(t)] 

+Bn/0E[X02(0] 

+B0tnE[x2n(t)] 

~4>an  (5.42f) 

The  parameters  in  equations  (5.19),  (5.25)  -  (5.27)  are  as  follows: 
For  the  case  fd(T)  =  \i£~& , 

/(p(0+y(*)/Vc)-/<exp{-{*(js(0+y(0/Vc)} 
-/(W0)+ft(W0);jj£ 

where 

i'-ic 

Similarly,  if  A0(T)  =  X0eVoT t  then 

4(W)  +  n0/Vc)  =  Aoexp{^(/?(f)  +  Y(0/Vc)} 

=a;/0(^w)+a;/1(«o)^ 


where 


1o  =  !7oC 


Also 
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A"n(/3(0  +  Y(t)  I VC)  =  A„[exp{r,;(/3(0  +  Y(()  /  Vc)}  - 1] 
=  *;,o(/5(<)Mnl(/J(f))^ 


where 


In  =  VnC 


A;,,(/5('))  =  A„r,>^C). 


*0  —   K0(- 


Kn  -  Krf- 


T(()=T(I)/C 

v  =  vC. 
A  steady-state  solution  to  equations  (5.39a)  -  (5.39c)  for  a  constant  toxin  input 

lit)  =  t  would  satisfy 


0  =  lie-Z'P{l-ao-an)-<pan-Xn[e^-l)an 


0  =  T*  - r-  U^«n  "  ~ *~  l>o/to0 

where  we  are  using  functions  of  the  form 

X0{T)  =  l0e*»T 
An(T)  =  A„(e"»T-l) 

with 


(5.43a) 
(5.43b) 

(5.43c) 
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T\o  =  t\(££  =  &,  r\n  =  VnC'  Vn  =  VnC,  V0  =  V0C,  Kn  =  KnC,  K0  =  K0C,T    =  T  /  C . 

Simplification  of  (5.43a)  and  (5.43b)  results  in 

A0 


an=ne  ZP 


^nfe^-l|+0  + 


fie 


■a 


-i-l 


(5.44a) 


(5.44b) 


Finally  expression  (5.44a)  and  (5.44b)  can  be  substituted  into  equation  (5.43c)  to 
obtain  an  equation  for  p.  This  latter  equation  can  have  0, 1,  or  2  solutions.  If  it  has 
no  solution,  then  an  =  Og  =  0  and  the  organ  is  dead.  If  the  equation  has  two 
solutions,  then  the  smaller  of  the  two  possible  solutions  is  biologically  plausible 
since  the  smaller  solution  is  an  increasing  function  of  r,  whereas  the  larger 
solution  is  a  decreasing  function. 
Numerical  Examples 

The  following  parameter  values  are  used  in  the  numerical  examples  below. 
v0  =  0.05,  \i  =  0.5,  rj0  =  1,  {  =  0.5,  k0  =  1,  X0  =  0.08,  vn  =  0.1,  r/n  =  0.5,  Kn  =  1, 
Ah  =0.05,  0  =  0.1,  C  =  107. 

Figures  6-9  present  results  for  the  case  in  which  t(0  =  0  for  0  <  t  <  2  and  then 

r(/)  =  4  x  105  for  t  >  2.  The  moments  are  started  at  their  steady  state  values  for  no 

toxin  input.  Figure  6  presents  the  mean  number  of  old  cell  pairs,  C  x  (XoQ),  and 

,     nl/2 


the  standard  deviation  of  the  number  of  old  cell  pairs,  CxElX^  (t)        .  Note  that 

the  mean  number  of  old  cell  pairs  decreases  to  a  new  steady  state  value  below 
the  value  for  no  toxin;  the  standard  deviation  increases,  then  decreases,  then 
increases  again  to  a  new  steady  state  value  which  is  higher  than  the  value  for  no 
toxin.  Figure  7  presents  the  mean  number  of  young /new  cell  pairs,  C  x  an(t),  and 
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The  mean 


the  standard  deviation  of  young/new  cell  pairs    CxE[x„(f)] 

number  increases  to  a  steady  state  value  larger  than  that  for  no  toxin  input;  the 
standard  deviation  initially  decreases,  then  increases  and  finally  decreases  to  a 
new  steady  state  value  below  the  value  for  no  toxin.  Figure  8  presents  the  mean 
number  of  totally  differentiated  cell  pairs,  C  x  («o(0  +  anit)),  and  the  standard 
deviation  of  the  number  of  totally  differentiated  cell  pairs 

,-nl/2 


Cx 


xoz(0+2Exmxn(0 


+  E 


2(0] 


The  mean  number  decreases  to  a  new  steady  state  value  below  the  value  for  no 
toxin;  the  standard  deviation  increases  to  a  new  steady  state  value.  Figure  9 
presents  the  mean  amount  of  toxin,  C  x  pit),  and  the  standard  deviation 

r  2    il1/2 

CxE  T  (t)        .  Both  values  increase  to  a  steady  state  value. 

Figures  10-  13  present  results  for  a  large  input  of  toxin  for  a  short  time 
period:  approximately  a  bolus  input.  The  input  of  toxin  is  Tit)  =  40  x  105  for 
2  <  /  <  3  and  Tit)  =  0  otherwise.  Figure  10  presents  the  mean  number  of  totally 
differentiated  cell  pairs  and  the  standard  deviation  of  totally  differentiated  cell 
pairs.  The  mean  number  initially  decreases,  then  increases  to  the  steady  state 
value  with  no  toxin.  The  standard  deviation  initially  increases,  then  decreases  to 
the  steady  state  value  with  no  toxin.  Figure  11  displays  the  mean  number  of  old 
cell  pairs  and  the  standard  deviation  of  old  cell  pairs.  The  mean  number  initially 
decreases,  then  returns  to  the  steady  state  value  with  no  toxin.  The  standard 
deviation  initially  decreases,  then  increases,  then  returns  to  the  steady  state  value 
with  no  toxin.  Figure  12  displays  the  mean  and  standard  deviation  of  the  number 
of  young/new  cell  pairs.  The  mean  number  initially  decreases,  then  increases, 
after  which  it  returns  to  the  steady  state  value  with  no  toxin.  The  standard 
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deviation  exhibits  similar  behavior  but  on  a  different  scale.  Figure  13  presents  the 
mean  and  standard  deviation  of  the  amount  of  toxin. 

Figures  14-17  present  results  for  the  same  parameter  values  and  same  toxin 
input  as  Figures  10  -  13  except  that  <p  =  0.01  instead  of  0.1;  that  is,  the  mean  time 
until  a  young/new  cell  pair  becomes  an  old  cell  pair  is  100  =  1/0.01  rather  than 
10  =  1/0.1.  This  should  increase  the  number  of  young/nev.r  cell  pairs.  Figure  16 
shows  that  the  mean  number  of  young /new  cells  is  indeed  larger.  Figure  14 
presents  the  mean  and  standard  deviation  of  the  number  of  totally  differentiated 
cell  pairs.  Note  that  the  mean  initially  decreases,  then  increases  to  overshoot  the 
steady  state  value  with  no  toxin.  The  mean  then  decreases  to  the  steady  state 
value  with  no  toxin.  This  behavior  of  overshooting  the  steady  state  value  with  no 
toxin  has  been  observed  in  experimental  studies;  Portier  (1993). 
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APPENDIX  A 
Isham-Whittle  Theory 

In  this  Appendix  we  present  another  approach  to  approximating  the  non- 
linear expectations  in  equations  (5.3)  and  (5.4).  We  follow  the  example  of  Isham 
(1991)  and  Whittle  (1957)  by  (a)  writing  equations  for  second  moments  and  joint 
transforms,  and  (b)  "closing"  the  system  of  lower-moment  equations  by 
assuming  that  (D(t),  T(t))  is  bivariate  Gaussian;  we  call  this  approach  Gaussian 
closure  (GC).  The  plausibility  of  this  assumption  follows  from  Model  D.l  above. 
Asymptotic  perturbation  methods  in  section  5  will  show  that  a  bivariate 
Gaussian  process  is  of  natural  occurrence  if  the  system  of  cells  is  large,  as  is  true 
in  practice. 

To  proceed,  we  write  second-moment  equations: 

E[D2(t  +  A)|D(<),r(f)]  =  (D(f)  +  lfn(T(t))(C  -  D(())A  +  (D(/)-  l)2A(T(())D(»)A 

+D2 (()[l  -  (ft  (T(t))(C  -  D(t))A  -  X (T(( ))D(t) A)]  +  o( A). 

Simplification  and  limit-taking  provides 


dt 


D2  (0]  =  2E[D(0/i(T(0)(C  -  D(f ))]  -  2e[d2  (f)A  (T(t ))] 
+4n(T(0)(C-D(0)]  +  E[A(r(0)D(0] 


(A.l) 


E[r2(*  +  A)|D(f),T(0]  =  (T(0  +  l)2T(0A 


v   w      '      l  +  xT(t)         K  \      {  w  1+kT(0    J 


+  °(A) 
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from  which 


«-E 
dt 


T2{t)]  =  2z(t)E[T{t)]+  r{t)-2vE 


D(t)T2(t) 
l+xT(f) 


+  vE 


D(t)T(t) 

l+KT(t) 


(A.2) 


Next 


E[D(t  +  A)T(/  +  A)|D(0,r(0]  =  (D(t)  +  l)(T(0)/i(T(0)(C  -  D(/))A 

+(D(f)  -  l)(T(f))A  (T(f))D(0A  +  D(f)(T(f )  +  l)r(f)A  +  D(0(T(f )  - 1)  ^fiW)A 


+D 

which  leads  to 


(t)T(t)  1  -»{T{t))(C  -  D(0)  A  +  A  (T(O)D(OA  +  <t)A  +  ^  j^y 


) 


^E[D(Or(0]  =  E[r(0/i(r(f))(c-D(f))]-E[r(OA(T(f))D(0] 


<*/ 


+  T 


(0E[D(f)]-uE 


DM 


D(f)T(Q 
l+*T(f) 


(A.3) 


In  order  to  go  further,  it  is  necessary  to  parameterize  the  rate  parameters  and 
the  Michaelis-Menten  (M.-M.)  term.  We  set  as  before,  for  £,  jj.  >  0, 

p(T)  =  Ai0e-^/A(T)  =  VT?r-  (A.4) 

If  (D(r),  T(r))  has  a  bivariate  normal  distribution,  then  the  distribution  is 
determined  by  the  marginal  means,  the  marginal  variances  and  the  covariance. 
Other  moments  can  be  expressed  in  terms  of  these  qualities.  For  example 
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E[D2(t)T(t)]  =  E[D(t)fE[T(t)} 

+2E[D{t)]Cov[D{t),T{t)] 

+E[T(t)]Var[D(t)]  (A>5) 

Equations  (5.3)  -  (5.4),  (A.1)-(A.3)  can  be  evaluated  using  the  following 

approximate  numerical  procedure.  Put 

mD(t)  =  E[D(tj\,  mT(t)  =  E[T(t)],  mD2  (t)  =  E[D2(t)]  mj2  (t)  =  E[T2(t)},  and 

mDT(t)  =  E[D(t)T(t)]. 

mD{t  +  A)  =  mD(t)  +  A{n[mT{t))C  -  (n{mT(t))  +  X(mT{t)))mD{t)}       (5.3a) 

mT(t  +  D)  =  mT(t)+A\r{t)^v- ^-—-mDT{t)\  (5.4a) 

[  l+KmT{t)  J 

mD2{t  +  A)  =  mD2(t) 

+A{2/x(mT(0)CmD(0-2(/i(mT(0)  +  A(mr(0))mD2(0 

+/i(mr(f))C  +  (A(ifir(0)-/i(fffrW))'«D(0}  (A.la) 

mj2  (t  +  A)  =  mj2  (t) 

+A{2T(t)mT(t)+T(t)-2v         l         e[d(0T2(0] 

+ 7TmDT(0) 

l+ioiirW   Dru'  (A.2a) 

mDT{t  +  A)  =  mDT(t) 

+A{Ai(m7(0)Cmr(0-(A(mr(0)  +  /iK(0)>nDr(0 

l+KmT{t)    L  J  J  (A3a) 
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The  higher  order  moments  E[D(t)T2(t)}  and  E[D2{t)T(t)]  are  evaluated  using  the 
expressions  for  a  bivariate  normal  distribution  such  as  (A. 5). 
Comparison  of  the  Two  Systems  of  Moment  Equations. 

In  this  section  we  study  the  limit  as  C  — >  «>  of  Equations  (5.3)  -  (5.4)  and 
(A.l)  -  (A. 3)  using  the  scaling  of  the  transition  rates  (5.12a)  -  (5.12d).  The  limits  of 
the  equations  are  as  follows. 


jLp{t)=r*-valt)?lt) 
dtHK>  l+KP(t) 

+2BtdE[X(t)Y(t)] 


j-tE[Y2(t)]  =  At(t)  +  2BttE[Y2(t) 

+2BdtE[X(t)Y(t)] 

vKa{t)p(t) 


+2 ^^ 

[l+KP(tjf 


E\Y'(t) 


-2- 


vP(t) 

1+Kp(t) 


1- 


P(*)K 
l+K*P(t) 


(5.3b) 
(5.4b) 


(A.lb) 


E[X(t)Y(t)] 


(A.2b) 
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dt 


E[X(f)Y(0]  =  (B^  +  B„)E[X(f)Y(0] 

+BdtE[x2(tj\  +  BtdE[Y2(t)} 


v  k  a(t) 


1        J  fl  +  v, 


+E 


;[x(()y«)]{B(i-4ft+Ai]}  (A3a) 

The  fact  that  the  second  moment  equations  (A.la)  -  (A.3a)  contain  more  terms 
than  (5.25)  -  (5.27)  is  due  to  the  difference  between  a)  first  taking  the  limit  as  C  -4 
«>  for  the  joint  distribution  of  (D(0,  T(0)  and  then  finding  the  second  order 
moments  of  the  miting  bivariate  normal  distribution  and  b)  writing  the 
differential  equations  for  the  second  moments  and  then  taking  the  limit  as 
C  -> «».  In  the  latter  case,  it  appears  that  the  effects  of  higher  order  moments 
become  important. 

An  example  calculation  follows.  Rewriting  equation  (A. 3)  using  (5.8) 


dt    L^ 


E[(ca(0  +  C1/2X(0)(c^(f)  +  C1/2Y(0) 


=  E[n(cp(t)  +  C1/2Y(0)[Q3(0  +  C1/2  Y(f)][c(l  -  a(t))  -  C1/2X(f)} 
-e[a(cj3(0  +  C1/2V(0)[Cj3(f )  +  C1/2Y(0][ca(f)  +  C1/2X(f) 
+r(0E[Ca(f)  +  C1/2X(0] 


-a 


E  (Ca(0  +  C1/2X(0)  (Cp(t)  +  C1/2  Y(f )) 


1+KTl 


K(cp(t)  +  Cl/2Y(tj) 


(A.4) 


Next  scale  the  transition  rates  using  (5.12a)-  (5.12d)  and  divide  both  sides  of 
(A.4)  by  C.  Equation  (A.4)  becomes 
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ftE[ca(t  )P(t)  +  Cl/2{P(t)X(t)  +  a(t)Y(t)}  +  X(f)Y(0] 

=  E  ^(0(0)  +  Aii(/?(0)^]c/?(0(l -  a(0)  +  C1/2[(l - a(f ))Y(f ) - p(t)X(t)] - X(t)Y(t) 

Ao(«0)  +  A*i(«0)^}ca(0/*(0  +  C1/2[«(f  )Y(0  +  p(t)X(t)]  +  X(f)Y(r)] 
+T*(f)E[ca(f)  +  C1/2X(0] 


Ve 


no 


i+*/?W   (i  +  v>(o) 


2  Vc 


-i>  E 


Y(f) 


l+*«0     [UrnTftt))2^ 


Ca{t)2p(t)  +  2C1/2a(0/?(f  )X(f)  +  X(02/J(/)] 


VCa(02  Y(t)  +  2a{t)X{t)Y{t)  +  C-y2Y{t)X{tf 


(A.5) 


If  a(t)  and  /3(0  satisfy  (5.19)  then  the  equation  corresponding  to  terms  of  order 
C  in  (A.5)  is  satisfied.  The  equation  corresponding  to  terms  of  order  C1/2  in  (A. 2) 
is  satisfied  if 

e[x(o]  =  E[y«)] = j-t  £[x(0]  =  |  £[>'«)]  =  o. 

The  equation  corresponding  to  terms  of  order  C°  in  (A.5)  is  equation  (A.3a). 
Equations  corresponding  to  terms  of  other  orders  become  0  as  C  — >  <». 


45 


APPENDIX  B 
Stochastic  Fluid  (Brownian  Motion)  Model  of  Toxic  Level 

The  Markov  chain  model  of  (5.1)  can  be  made  more  plausible  and  flexible  by 
replacing  the  discrete  state  space  by  a  continuous  one.  It  is  possible  to  allow 
[T(t)},  hereafter  called  toxicity,  to  move  according  to  a  Brownian  motion  or 
Wiener  process  with  drift;  both  drift  and  diffusion  coefficient  can  depend  on  the 
current  number  of  differentiated  cells  and  the  toxicity. 

Derivation  of  the  replacement  of  the  forward  equation  in  transform  form, 
(5.7),  by  the  consequence  of  the  above  T(t)  representation  gives  (5.7)  with  the  last 
two  lines  replaced  by 


6fE 


1  9 

2  l 


BMi>etnt)5{p[t)iT{t)y 
e^D^+6^t'a2{D(t)>T(t)) 


(B.l) 


here  8(D(t),T(t))  and  c2(D(t),  T(t))  are  the  Brownian /Wiener  state-dependent 
drift  and  diffusion.  In  turn,  these  must  be  scaled: 


5(D(t)J(t))  =  C8*Ut)  +  ^At)  +  ^ 


a2(D(0,nO)-Co2^a(0  +  ^«0+^ 


(B.2a) 


(B.2b) 


These  functions  can  then  be  Taylor-series  expanded  to  provide 
8{D(t)J(t))  =  C?(a(t),ftt^ 

o2(D(t)J(t))  =  CcS\a(t)/p(t))  +  ^(rt{a(t\ 

Now  apply  (5.9)  and  the  new  version  of  (5.11)  after  scaling  as  in  (5.12);  the 
change  occurs  in  the  last  two  lines  of  (5.11)  as  follows 
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^E[e9'xW+«'r(|)(ci,(«(0./J(0)  +  Vc(«1,1(a(/)./J(0)X(»)+51*2(«(»),/J(0)y(0))] 

+i^E[e^x(,)+e,y«)(Cc72-(a((W(/))+vc(ff?1-(a(0,/3«))x(() 
+og(o(0^(0)xW+oS(o(«).«OMO)) 

-Vce^-5*(a(/),«0)+e/«,,,(a(0.«0)#  +  ^2(tt(').«0)^-]  (B.4) 

+Ie,2pa2,(a((),«0)  +  O(l/Vc), 

where  <5u  and  5^2  are  partial  derivatives  of  &  with  respect  to  a  and  /?,  at  oit)  and 
0(0. 

From  (5.18)  with  the  coefficient  of  fy  changed  as  above,  (B.4),  we  obtain  the 
replacement  for  mean  scaled  toxicity, 

and  from  (B.4)  and  (5.20)  we  now  find 

Bdl(t)  =  5'n(a(t),m).    Bu(t)  =  5l2(a(t),P(t)).  (B.5) 

Example. 


(B.6) 


(B.7) 


of  and  of  are  (here,  but  not  necessarily)  constants  that  permit  adjustment  of 
variances  of  toxic  agent  input  I  of )  and  output  (of)-  Scaling  as  before, 
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8(D(t)T(t))   c.m    ^(o(MO+Vcx(0)(c/3(Q,Vcy(0)) 
5(d(0,t(0)-Ct  «— 1  +  ^(o+Y(0/Vc) 

((^)+x(0/Vc)(W)+y(Q/Vc)) 
=  C(T  ^  (°)" i+/(W)+yW/Vc) 


02(0(0^(0)  =  c 


Bc**(a(0+x(0/Ve,Y(0/V£)  (B.8) 

t/  (n((«(f ) + x(0  /  Vc)(/?(0 + Y(f )  /  Vc ))" 


0?T*{t)+0% 


J 


=  cc?(a{t) + x(0  /  Vc,/?(0 + y(0  /  Vc)  (B  9) 

From  these  we  get 

^-  =  -v(t)      PVal'£MtW))  =  Bil(t).  (B.10) 

<?a  1  -H  KT  /3(f) 

af/(0+^^(0-^^-^(«(0,«0)  =  AW-  (B.12) 

Here  ^o(/^(0)  =  K  0(0  an<^  vl(0  =  ^  m  (5-21);  the  agreement  with  the  previous 
model  is  apparent  when  of  =  of  =  1.  Note  in  particular  that  the  model  (5.1) 
represents  toxic  chemical  input  as  a  (time-dependent)  Poisson  process.  To 
represent  a  deterministic  inflow  simply  put  of  =  0,  while  to  represent  an  extra- 
Poisson  variability  put  of  >  1.  Variability  in  the  output/removal  of  toxicity  can 
be  similarly  modeled  by  adjusting  of. 

Figures  IB,  2B,  3B  present  the  standard  deviation  of  the  number  of  totally 
differentiated  cell  pairs  for  of  =1,  5,  and  20.  Figure  IB  has  a  constant  input  of 
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toxin.  Figures  2B  and  3B  have  a  pulse  of  toxin  but  different  values  of  4>.  As  of 
increases,  the  standard  deviation  of  the  number  of  totally  differentiated  cell  pairs 
also  increases. 
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